↜ Back to index Introduction to Numerical Analysis 1

Today we focus on programming that involves loops: executing the same code multiple times.

do loops in Fortran

In Fortran, the basic loop construct is the do loop:

! prints numbers 1, 2, ..., 5
implicit none
integer i

do i = 1, 5
    print *, i
enddo
end

outputs

           1
           2
           3
           4
           5

Fortran will run the code between do and enddo statements for every value in the provided range (1, 5 stands for 1, 2, 3, 4, 5) assigned to variable i. The above code is equivalent to

! prints numbers 1, 2, ..., 5
implicit none
integer i

! unrolled do loop
print *, 1
print *, 2
print *, 3
print *, 4
print *, 5

end

You can add a third number in the range that specifies that some numbers should be skipped.

! prints numbers 1, 3, 5, 7, 9
implicit none
integer i

!             👇 step size (default 1)
do i = 1, 10, 2
    print *, i
enddo
end

outputs

           1
           3
           5
           7
           9

Note that 10 does not appear.

The loop variable (i in the above example) should be an integer. It is tempting to use a real variable at times, but this gives a warning when compiling.

implicit none
real x

do x = -1,1,0.5
    print *, x
enddo

end

$ gfortran code/do-loop-real.f90
code/do-loop-real.f90:4:3:

    4 | do x = -1,1,0.5
      |   1
Warning: Deleted feature: Loop variable at (1) must be integer
code/do-loop-real.f90:4:12:

    4 | do x = -1,1,0.5
      |            1
Warning: Deleted feature: Step expression in DO loop at (1) must be integer

The code still works if you run it, but it is probably better to use an integer variable for looping and compute the real value from it. You can also avoid weird issues caused by rounding of real numbers. The following code produces the same output, and Fortran is happy.

implicit none
integer i
real x

do i = 0, 4
    x = -1 + 0.5 * i
    print *, x
enddo

end

Loops can be nested. This is useful when we need to run a code for all pairs (triples, etc.) of values. It is a good idea to indent the code in each inner loop to make it easier to read.

! Prints the coordinates of the vertices of a unit cube
implicit none
integer x, y, z

do x = 0, 1
    do y = 0, 1
        do z = 0, 1
            print *, x, y, z
        enddo
    enddo
enddo
end

           0           0           0
           0           0           1
           0           1           0
           0           1           1
           1           0           0
           1           0           1
           1           1           0
           1           1           1

! Prints the truth table of logical OR and logical AND
! for all possible inputs.
implicit none
integer i, j, or, and

print *, '          a |           b |      a OR b |     a AND b'
print *, '------------------------------------------------------'
do i = 0, 1
    do j = 0, 1
        ! logical or: at least one is 1
        if (i + j >= 1) then
            or = 1
        else 
            or = 0
        endif
        ! logical and: both are 1
        if (i + j >= 2) then
            and = 1
        else 
            and = 0
        endif
        print *, i, "|", j, "|", or, "|", and
    enddo
enddo
end

           a |           b |      a OR b |     a AND b
 ------------------------------------------------------
           0 |           0 |           0 |           0
           0 |           1 |           1 |           0
           1 |           0 |           1 |           0
           1 |           1 |           1 |           1

! Prints a nice multiplication table.
! Also shows how to format output.
implicit none
integer a, b

do a = 1, 9
    do b = 1, 9
        !         👇I2: integer of width 2 characters
        !            👇A: a string (we use it to print the comma ',')
        !                👇advance='no': do not insert a new line
        write(*,'(I2,A)',advance='no') a * b, ','
    enddo
    ! new line
    write(*,*) '' 
enddo
end

 1, 2, 3, 4, 5, 6, 7, 8, 9, 
 2, 4, 6, 8,10,12,14,16,18, 
 3, 6, 9,12,15,18,21,24,27, 
 4, 8,12,16,20,24,28,32,36, 
 5,10,15,20,25,30,35,40,45, 
 6,12,18,24,30,36,42,48,54, 
 7,14,21,28,35,42,49,56,63, 
 8,16,24,32,40,48,56,64,72, 
 9,18,27,36,45,54,63,72,81, 

It is possible to stop a loop early using the exit statement. This following pattern is common when searching for the first value that satisfies a certain condition.

! Computes the square root of an integer by trying all possible values.
! (Not a great algorithm... use built-in function sqrt :) )
implicit none
integer i, n, c

print*, 'Enter an integer'
read*, n

c = 0   ! initialize the flag
do i = 0, n
    if (i * i == n) then
        print*,'Square root of', n, 'is', i
        ! No need to continue checking, we stop the loop.
        ! We keep a **flag** to avoid printing a message after the loop.
        c = 1
        exit
    endif
enddo

if (c == 0) then
    ! value 0 indicates that we haven't found a square root
    print*, n, 'does not have an integer square root'
endif

end

Common do loop algorithm

Loops are great for doing calculations with a sequence of numbers.

• sum a_1 + a_2 + \cdots + a_n
• product a_1 \cdot a_2 \cdots a_n
• maximum \max(a_1, a_2, \ldots, a_n), minimum \min(\ldots)
• find nth value s_n of a recurrence formula s_i = f(s_{i-1}, a_i), i = 1, \ldots, n for given s_0.

As a concrete example, let us consider the finite sequence a_i = i(5 - i), \qquad i =1, \ldots, 10

Let us start with the sum. How do we express the sum a_1 + a_2 + \cdots + a_n as a sequence of commands with a common pattern? Since the order of additions does not matter, we usually perform them from left to right:

(\cdots (((a_1 + a_2) + a_3) + a_4) \cdots +a_n)

So we are performing a sequence of additions \begin{aligned} s_2 &= a_1 + a_2\\ s_3 &= s_2 + a_3\\ s_4 &= s_3 + a_4\\ &\vdots\\ s_n &= s_{n-1} + a_n \end{aligned}

And s_n is the final answer. We can see a nice pattern in the sequence of computations. The first one s_2 = a_1 + a_2 is a bit special, but the rest is just s_i = s_{i-1} + a_i.

When writing this as a sequence of commands for a computer, we want to first think about what values we need to remember so that we can declare necessary variables.

• a_i are given by a formula in our case, or come from an array, or some other data source (read from the keyboard, etc.). Depends on what we are computing a sum of.

• s_i: Here we note that we need each s_i only in the very next computation. It would be wasteful to keep them all in variables. In fact, we need only one variable to keep the latest computed s_i, which we then overwrite. We could change the above computation by this “pseudocode”:

  \begin{aligned} s &\leftarrow a_1 + a_2\\ s &\leftarrow s + a_3\\ s &\leftarrow s + a_4\\ &\vdots\\ s &\leftarrow s + a_n \end{aligned}

  Here \leftarrow is the command “assign the value on the right to the variable on the left”. This is the variable = value Fortran statement. The meaning is very different from a = b used in math.

With this we can write the code that computes the sum for a_i = i(5-i), i = 1, \ldots, 10 in our example.

! Sum: version 1
implicit none
integer s

s = 1 * (5 - 1) + 2 * (5 - 2)
s = s + 3 * (5 - 3)
s = s + 4 * (5 - 4)
s = s + 5 * (5 - 5)
s = s + 6 * (5 - 6)
s = s + 7 * (5 - 7)
s = s + 8 * (5 - 8)
s = s + 9 * (5 - 9)
s = s + 10 * (5 - 10)

print *, s
end

This works, but it is rather annoying. What if we want to compute the sum of million of terms, or the number of terms depends on the size of input data? Of course we should use a do loop. We introduce a loop variable, say i, and put the common statements in a loop.

! Sum: version 2
implicit none
integer s, i

s = 1 * (5 - 1) + 2 * (5 - 2)
do i = 3, 10
    s = s + i * (5 - i)
enddo

print *, s
end

Much better. The computer does exactly the same sequence of commands and in the same order as before, but it is much more readable and we need to write the formula for a_i only 3 times. Can we do even better?

Yes. We can rewrite the step s \leftarrow a_1 + a_2 as \begin{aligned} s &\leftarrow 0\\ s &\leftarrow s + a_1\\ s &\leftarrow s + a_2 \end{aligned}

This simplifies the first step of the algorithm to s \leftarrow 0 and the two following steps can be then included in the loop.

! Sum: version 3
implicit none
integer s, i

s = 0
!     👇 starting from 1 now
do i = 1, 10
    s = s + i * (5 - i)
enddo

print *, s
end

What about the maximum \max(a_1, a_2, \ldots, a_n)?

To use the same algorithm as before, we observe that we can write it as \max(\ldots\max(\max(a_1, a_2), a_3), \ldots, a_n)

so we can compute the maximum as \begin{aligned} m &\leftarrow \max(a_1, a_2)\\ m &\leftarrow \max(m, a_3)\\ m &\leftarrow \max(m, a_4)\\ &\vdots\\ m &\leftarrow \max(m, a_n)\\ \end{aligned} The first step can be rewritten as \begin{aligned} m &\leftarrow a_1\\ m &\leftarrow \max(m, a_2)\\ \end{aligned}

We could also try to rewrite it as \begin{aligned} m &\leftarrow -\infty\\ m &\leftarrow \max(m, a_1)\\ m &\leftarrow \max(m, a_2) \end{aligned}

However, -\infty is a not a valid value of an integer in Fortran. We could replace it by a “sufficiently” small number instead or even the smallest possible integer -2^{31}. This depends on the sequence.

Finally, it is left to implement m \leftarrow \max(m, a_i) in Fortran. Here we can use the if condition

if (ai > m) then
    m = ai
endif

or the built-in max function:

m = max(m, ai)

Here is the complete code:

! Find the maximum of a simple sequence.
implicit none
integer m, ai, i

! set m to a_1
m = 1 * (5 - 1)
!      👇 starting from 2
do i = 2, 10
    ai = i * (5 - i)
    if (ai > m) then
        m = ai
    endif
enddo

print *, m
end

Sometimes we want to know which element is largest, that is, the value of i such that a_i = \max(a_1, \ldots, a_n). This is sometimes called arg max. There can be multiple such i, in which case we need to decide whether to find all of them, or only the first one, …

We can modify the above algorithm to keep track of for which i we update m. This will return the smallest i so that a_i is equal to the maximum.

! Find the index of the largest element: argmax
implicit none
integer m, ai, i, im

! set m to a_1
m = 1 * (5 - 1)
! position of the largest element is initialized to 1
im = 1
!      👇 starting from 2
do i = 2, 10
    ai = i * (5 - i)
    if (ai > m) then
        m = ai
        ! remember i when m is updated
        im = i
    endif
enddo

print *, im
end

Summary

The algorithm in this section can be used whenever the computation can be written as an initialization step s \leftarrow value followed by a sequence of steps s \leftarrow f(s, i).

Next time we will see that some numerical methods for solving ordinary differential equations are in this form.

Short report problems

$ ./a.exe
 Enter a positive integer
12
           6
           3
          10
           5
          16
           8
           4
           2
           1
           4

$ ./a.exe
 Enter a positive integer
12
          16

$ ./a.exe
 Enter a positive integer
8
           4

Other exercises

$ ./a.exe
 Enter an integer
8
 Not a prime

$ ./a.exe
 Enter an integer
1
 Not a prime

$ ./a.exe
 Enter an integer
2
 Prime number

$ ./a.exe
 Enter an integer
93
 Prime number